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Abstract 

Despite recent advances in the field of Oncoimmunology, the success potential of immunomodulatory 
therapies against cancer remains to be elucidated. One of the reasons is the lack of understanding 
on the complex interplay between tumor growth dynamics and the associated immune system responses. 
Towards this goal, we consider a mathematical model of vascularized tumor growth and the corresponding 
effector cell recruitment dynamics. Bifurcation analysis allows for the exploration of model’s dynamic 
behavior and the determination of these parameter regimes that result in immune-mediated tumor control. 
Here, we focus on a particular tumor evasion regime that involves tumor and effector cell concentration 
oscillations of slowly increasing and decreasing amplitude, respectively. Considering a temporal multiscale 
analysis, we derive an analytically tractable mapping of model solutions onto a weakly negatively damped 
harmonic oscillator. Based on our analysis, we propose a theory-driven intervention strategy involving 
immunostimulating and immunosuppressive phases to induce long-term tumor control. 


Introduction 

The immune system is widely recognized for its capacity to detect and destroy cancer cells, as well as 
to prevent tumor recurrence maintaining an immunological memory 1 . Indeed, every known innate 
and adaptive immune effector mechanism has been reported that participates in tumor recognition and 
rejection |2j. However, experimental observations support that tumorigenic processes themselves can 
promote immunosuppression or immune tolerant states that facilitates neoplastic growth and progression 
J3}[4| . Cancer cells employ diverse strategies to inhibit or block immune responses, including tumor- 
induced impairment of antigen presentation, secretion of immunosuppressive cytokines or expression of 
surface molecules, as well as production of diverse pro-apopoptic factors 5]. Nevertheless, there are 
clinical and preclinical evidences supporting that activation of the innate antitumor immunity can result 
in tumor regression and provide therapeutic benefits |6 . 

The main goal of oncoimmunology is to strengthen the immune system’s innate ability to combat 
and kill cancer cells by enhancing the effectiveness of the immune responses. Among the different im- 
munotherapeutic techniques are checkpoint inhibitors, immune response modifiers (cytokines), mono¬ 
clonal antibodies and vaccines [7]. Passive and active immunotherapy has been successfully applied to 
the treatment of a wide variety of human cancers [8] and holds promise of a lifelong cure [9 . However, 
tumor-induced immunosuppression still represents a major obstacle to effective cell-mediated immunity 
and immunotherapy [dj[5]. Accordingly, more insights into the main mechanisms associated with im¬ 
mune responses based on tumor specific features are required to obtain successful therapeutic outcomes 
with immunomodulatory strategies. Despite years of research devoted to understand the underlying 
mechanisms of immune-tumor interactions, there are still many unanswered questions. In particular, 
those related with the impact of tumor-associated vascularization on immune responses, as well as de¬ 
termination of optimal and effective therapeutic protocols in cancer immunotherapy, are far from being 
completely elucidated. 
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Mathematical oncology is a valuable descriptive and predictive analytic framework to address such 
open questions. Continuum Mechanics concepts have been widely considered to investigate tumor growth 
and therapy implications 10 14 . In addition, several mathematical models of tumor growth, where some 


forms of the immune dynamics are often included, have been also extensively studied in the last years 
15-24 . Clinical data evidence that cancer cells can survive in a undetectable dormant state for extended 


periods of time 25 
20 22 20:27 


which has been also predicted by several models of tumor-immune cell interactions 19 


However, the neoplasm develops diverse strategies to circumvent the anti-tumor action of the 
28p9 . In particular, this equilibrium state can be disrupted by different events affecting 


immune system 

the immune system which could result in tumor regrowth 29 


system have been observed both in its healthy state and pathological situations 


Sustained oscillations by the immune 
27p0pl 


Therefore, the 


presence of an immune component in mathematical modeling has been described crucial for reproducing 
clinically observed phenomena such as tumor dormancy, oscillations in tumor size and spontaneous tumor 
regression 27 32 -36 . Among the several reviews on the subject are those covering mathematical models of 
tumor growth mainly focused on the cancer and immune system interactions 37 -45 . The mathematical 
modeling of the entire immuno-oncology dynamics is an enormously difficult and complex task. In 
consequence, models describing interactions between growing tumors and immune dynamics should focus 
on the crucial factors that are known to allow tumor escape from immuno-surveillance. Tumor-induced 
angiogenesis is a crucial mechanism for cancer survival and proliferation, allowing a continuous supply 
of oxygen and nutrients needed for tumor growth and progression 46-48 . However, the effectiveness of 


antitumor immune responses is associated with the functional levels of the tumor blood vessels, which 
allow a wider range of effector cell types to penetrate the tumor bulk and further exterminate cancer 
cells [49|50 . These opposing effects demand for a mathematical model of vascularized tumor growth that 
allows to explore the therapeutic potential of immunomodulatory interventions when innate immune 
responses are insufficient for long-term tumor control. 

The work herein reported intends to yield new insights in the potential of immunomodulatory inter¬ 
ventions for cancer therapy. To this end, we propose a tumor-effector cell model based on well-known 
biological assumptions that combines a model of radially symmetric tumor growth with an immune cell 
recruitment model [20 p2pl| |52| . The main feature of our model is the modeling of the interplay between 
functional tumor-associated vasculature and effector cell dynamics as described in 53]. Model results 
predict that, depending on the functional tumor vascularization degree and effector cell recruitment 
rate, long-term tumor control cannot be always reached. We particularly focus in such situations where 
tumors escape immuno-surveillance to suggest an optimized theory-driven therapeutic strategy against 
tumor growth. A temporal multiscale approach is then implemented to describe the tumor-immune sys¬ 
tem interactions, where an analytically tractable approximation of the cancer-effector cell dynamics is 
derived. We find that an efficient modulation of the immunostimulating and immunosuppressive phases 
could induce long-term tumor control. 


Mathematical model description 


The present model describes the interactions between growing tumors and induced immune system dy¬ 
namics. More precisely, the proposed tumor-effector cell model is a combination of a radial tumor 


growth model and an effector cell recruitment model originally proposed in 51 and 32 respectively, see 
also 20 52 54 55 . The main feature is the low dimensional modeling of the complex interplay between 
tumor-associated functional vasculature and immune recruitment dynamics as described in |53| . This 
model can be interpreted as the temporal evolution of the average tumor radius, since radially symmetric 
growth is not a realistic behaviour. The system variables are the average tumor radius R(t) and effector 
cell concentration E(t) in the tumor vicinity. 
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Radial tumor growth, R(t) 

The temporal evolution of the average tumor radius is considered, where for simplicity invasive and 
diffusive tumor properties are not taken into account. The tumor is modeled as an incompressible fluid 
flowing through a porous medium, where tissue elasticity is simplified. The tumor-host interface is 
assumed to be sharp and cell-to-cell adhesive forces are modeled as a surface tension at that interface. 
The tumor expands as a mass whose growth is governed by a balance between cell birth (mitosis) and 
death (apoptosis and necrosis). The mitotic rate within the tumor is assumed to be linearly dependent 
on the nutrient concentration (oxygen, glucose, etc.) and is characterized by its maximal value A m 
at the tumor-host interface. The death rate is uniform within the tumor and constant in time. 
Moreover, we assume that the death rate A^ reflects the lump effect of apoptotic/necrotic processes and 
any other cell death inducing factor. The concentration of nutrients (e.g. oxygen or glucose) obeys a 
reaction-diffusion equation in the tumor volume, where nutrients are supplied from the tumor-associated 
functional vasculature and consumed by the tumor cells at a uniform consumption rate. 

To gain insight on the impact of vascularization in tumor growth and immune responses, we assume 
that the non-negative and dimensionless parameter B , where 0 < B < 1, represents the net effect of 
functional tumor-associated vasculature on the tumor radius evolution. In the limit of avascular tumor 
growth B = 0, such tumor-effector interactions take place only at the tumor surface. At the other 
extreme, for B — 1 effector cells can potentially interact with any cancer cell within the tumor bulk. 
Moreover, we consider an intrinsic length scale Lp representing the average length of nutrient gradient, 
i.e. supply, diffusion and consumption. 

The efficacy of immune killing depends on the ability of effector cells to penetrate the tumor bulk via 
the functional tumor-associated blood vessels. With improved vascularization, the effectors kill tumor 
cells not only on the surface of the tumor, but also further inside [56 . We consider this process through a 
phenomenological scaling function /(!?,!?) = r b-i +1 G [0,1], for B G [0,1), that models the penetration 
of effector cells in the tumor parenchyma through the existing functional vasculature. This function 
modulates the term related to tumor-effector cell interactions, such as killing of tumor cells due to 
effectors represented by a rate c. Such scaling functions have been also considered in the classical von 
Bertalanffy approach, and more recently in allometric models |57| . 

Taking these factors into account, we deduce the tumor radius R(t) dynamic under the assumption 
of radial symmetry according to the following ODE equation: 

f = 1(A mB - , a)R + A„(l - B)L„ ( tanh( ; /Lg) - i §) - cERH^m (1) 

net vascular avascular death due to 

tumor growth tumor growth effector cells 

where the time coordinate t has been omitted for notational simplicity, and A m-, A^, Le> and c are 
non-negative constants. 


Effector cell concentration, E(t) 

We assume that effectors are recruited at a rate r depending on tumor cells following the Michaelis-Menten 
kinetics, where K is a positive constant denoting the concentration at which the immune recruitment is 
half-maximal. Effector cells die at a rate d 0 , and become inactivated a rate di due to their antitumor 
activity. In particular, the inactivation of effectors by tumor cells is modeled through the function f(R, B) 
described above. As functional vascularization increases, effectors can kill cancer cells throughout the 
tumor bulk and tumor-immune cell interactions increase resulting in more inactivated immune cells. 
Moreover, innate immunity or base immuno-surveillance is represented as a minimum presence of active 
effector cells at any time given by an innate immune growth rate a, even in the absence of tumor cells. 
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The resulting ODE equation for the effector cell concentration E dynamic is given by: 


E - diER 3 f(R, B) - d 0 E + 



( 2 ) 


immune effectors effectors innate 

recruitment rate inactivation rate death rate immunity rate 

where the time coordinate t has been omitted for sake of simplicity in the notation, and r, K , do, e?i 
and a are non-negative constants. 


Model parameterization 

Under the small tumor radius assumption, the very early tumor growth is always of exponential nature 
and does not depend on the vascularization effects, i.e. parameter B , 


52 54 58 59]. Accordingly, we 
assume that this initial growth depends exclusively on the net proliferation rate X p = (Am — A^) in the 
absence of adaptive antitumor immune responses at those stages of growth, see the first term of Eq. 0 - 
Thus, considering experimental estimates of the growth rate at early phases of spheroids tumor growth 
for the mouse colon carcinoma cell line CT26 as A p ~ 1.20 day -1 
value Am = 1/18 ft. = 1.34day -1 


for CT26 murine cells 


53 


60 


61 


53 


60 


and the physiological plausible 
we have that A^ ~ 0.14day -1 . The 


characteristic nutrient diffusion length has been experimentally estimated that ranges between 0.2 and 
0.3 mm 53,62-64 . The effector cell characteristic concentration is at the order of magnitude 10 5 cells. 


The latter estimate is justified since the characteristic length scale of the system is at the order of 1.0 mm, 
and given that cells are commonly assumed with a diameter between 10/im and 20/rm [611, then for a 


volume of 1mm 3 the concentration is at the order of magnitude 10 5 cells. Moreover, we consider c = 0.03 
cells -1 days -1 as measured from murine CT26 tumor growth experiments in |53| , see also [20p2]|55] . The 
remaining parameter values do = 0.37 day -1 , d\ = 0.01 mm -3 days -1 , K = 2.72 mm 3 and a = 0.13 x 10 5 
cells days -1 are considered from previously reported experimental data |20 ,[32][55jj65j|66|, and properly 
rescaled to the magnitudes and units considered in our model. Through a parameter sensitivity analysis, 
we explore the effects on tumor growth of the effector cell recruitment rate r and functional vascularization 
degree B for different choices of the initial tumor radius Rq and concentration of effector cells Eq- 

For convenience of the reader, we summarize in Tab. |T] the parameter values used in numerical simu¬ 
lations of the tumor-immune dynamics considering the effect of tumor-associated vasculature. 


Description Parameter Value Units Sources 


Tumor mitotic rate 

Am 

1.34 

days -1 

53 

60 

61 


Tumor death rate 

Aa 

0.14 

days -1 

53 

60 

61 


Characteristic nutrient diffusion length 

Ld 

[0.2 - 0.3] 

mm 

53 

62 

04 


Tumor cell kill by effectors 

C 

0.03 

cells -1 days -1 

20 

32 

53 

55 

Tumor volumen where r is half-maximal 

K 

2.72 

mm 3 

20 

32 

55 

66 

Effectors inactivation rate by tumor cells 

di 

0.01 

mm -3 days -1 

20 

32 

55 

66 

Effectors death rate 

do 

0.37 

days -1 

20 

32 

65 

66 

Innate immune growth rate 

a 

0.13 x 10 5 

cells days -1 

20 

32 

55 

66 


Table 1: Model parameters. The effects of model parameters B and r, i.e. functional tumor-associated 
vasculature and effector cell recruitment rate respectively, are investigated by a model sensitivity analysis. 


Dynamical analysis of the model 

In this section, we analyze the model’s behaviour with respect to two parameters, namely the effector 
cell recruitment rate r and functional tumor-associated vasculature B. 
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Fixed point analysis 

The first step towards analyzing the model dynamics is the identification of the fixed points along with 
their stability classification. Fig. [lj A-D) depicts the phase portraits of the system of equations <[l])-([2j) for 
different values of the model parameter B , while keeping r constant and equal to 0.57 days -1 . The black 
curves represent the nullclines, i.e. curves along which dR/dt = 0 and dE/dt = 0, and the colored curves 
the system trajectories for different initial conditions. The fixed points of the system are located at the 
intersection of the nullclines. In each case, we identify the existence of two fixed points corresponding to 
a low tumor radius (L = (Rl, El)) and a high tumor radius (H = (Rh, Eh))- 


Eo = 20 x io 5 cells 


Eo = 10 x 10 5 cells 


A Vasculature (B = 0.50) B Vasculature (B = 0.65) C Vasculature (B = 0.80) 



D Vasculature (B = 0.95) 




Figure 1: Phase portraits and long-term characteristic dynamics. (A-D) Phase portraits of 
the tumor radius R versus the concentration of effector cells E for different functional levels of the 
tumor-associated vasculature B. The colored lines depict trajectories starting with the initial conditions 
IC\ = (Rq, Eo) = (2.0 mm, 20 x 10 5 cells) (green) and IC 2 = (Rq,E 0 ) = (2.0 mm, 10 x 10 5 cells) (purple), 
respectively. The nullclines (zero-growth isoclines) of the dynamical system (black lines) are also plotted. 
(E-H) Temporal evolution of R (red lines) and E (blue lines) which correspond to the trajectories in 
panels (A-D) starting at the initial condition IC\. (I-L) Temporal evolution of R (red lines) and E (blue 
lines) which correspond to the trajectories in panels (A-D) starting at the initial condition /C 2 . The 
immune recruitment rate is r = 0.57 days -1 and the remaining parameters are as in Tab. [ll 


Fig. HE -L) shows the evolution of the tumor radius and effector cells of the corresponding trajectories 
on the phase plain, presented in Fig. B A -D). A quick glance at Fig. [I] reveals that the system trajectories 
initiated near the L fixed point follow an oscillatory behavior with a slowly varying amplitude that can 
be either increasing Fig. 0F ,G) or decreasing Fig. [ljE,H). This behavior indicates that, depending on 
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the model parameter values, the L fixed point can be an attractor or a repellor. On the other hand, the 
trajectories initiated away from this point follow an exponential behavior which results in a boundless 
tumor growth while the amount of the effector cells fades out (see Fig. 0I-L))- 

To gain insight in the system behavior near the fixed points we construct the bifurcation diagrams with 
respect to the model parameter r for different values of B (see Fig. [2jA-C)). The bifurcation diagrams 
were calculated by performing an arc-length continuation method. This method is a special case of 
numerical fixed-point continuation methods that ensures the continuation of solution branches at turning 
points 


67 68 


For sufficiently small values of r, we obtain that there are no fixed points, which consequently implies 
that the tumor grows indefinitely. However, for a critical tumor radius r cr , a homoclinic bifurcation 
occurs 69 70 giving rise to two states: a lower branch which corresponds to the L fixed point, and an 
upper branch which corresponds to the H fixed point. The local stability analysis shows that the H 
fixed point is a saddle point, and therefore unstable. On the other hand, the L fixed point is a spiral 
point which, depending on the value of the parameter B , can be unstable (spiral source) or stable (spiral 
sink) m 72]. Fig.[2](D,E) shows the local stability analysis of L and H for B = 0.60, which corresponds 
to the bifurcation diagram in Fig. [2](B) . 


Vasculature (B = 0.40) 


Vasculature (B = 0.60) 


Vasculature (B = 0.95) 





Recruitment rate r (days- 1 ) 


Recruitment rate r (days- 1 ) 


Recruitment rate r (days- 1 ) 



Recruitment rate r (days- 1 ) 


xio - 3 E 



Recruitment rate r (days- 1 ) 


Figure 2: Bifurcation diagrams with respect to r and local stability analysis. (A-C) One- 
parameter bifurcation diagrams with respect to the effector cell recruitment rate r for different values of 
the functional tumor-associated vasculature B (0.40, 0.60 and 0.95, respectively). The upper branches 
correspond to the saddle point H, whereas the lower branches to the spiral point L. Solid lines depict 
stable fixed points and dotted lines the unstable fixed points. (D) Eigenvalues of the Jacobian estimated 
at the saddle point H with respect to the parameter r for B = 0.60. (E) Real part of the eigenvalues of 
the Jacobian estimated at the spiral point L with respect to the parameter r for B = 0.60. 


In the case of L being a spiral sink, a system trajectory initiated inside the homoclinic orbit, i.e. the 
closed orbit which starts and returns to the saddle point in Fig. [3]( A), will follow regular oscillations with 
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a slowly decreasing amplitude until it reaches the fixed point. This implies that the tumor radius will stay 
in a control bounded state. The homoclinic orbit defines the basin of attraction for the spiral sink (i!).70 


Thus, any trajectory initiated outside the homoclinic orbit will result in an uncontrollable tumor growth 
while the concentration of effector cells fades out. On the other hand, if L is a spiral source, any 
initialization of the system close to L will produce regular oscillations with a slowly increasing amplitude. 
This behavior persists until the trajectory of the system reaches the unstable manifold of the saddle point 
H which drives the system towards an exponentially uncontrollable tumor evolution (see Fig. |3j(B)). 


A 


40 



1 2 3 4 5 


Tumor Radius R (mm) 


B 



Tumor Radius R (mm) 


Figure 3: Classification of the system trajectories. (A) Homoclinic orbit of the saddle point H 
(star) when the spiral point L (dot) is stable. (B) The unstable manifold of the saddle point H (star) 
when the spiral point L (dot) is unstable 


Interestingly, the time required for reaching the unstable manifold is significantly large. This system’s 
behavior can be explained by performing a local stability analysis at the spiral point L. Fig. BA.B) 
illustrates the stability of L when r is kept fixed and B is allowed to vary. The eigenvalues of the 
Jacobian estimated at L are A = jj ± i/3. When /i < 0, L is a spiral sink, while for fi > 0, L is a spiral 
source. Fig. |4](C,D) shows how /i changes with respect to parameter B. For all values of the parameters 
considered, we find that when L is a spiral source the real part of the eigenvalues is /i < 1 with a 
maximum value of Umax — 5.5 x ICC 3 . Therefore, by setting fi = e <C 1, we can identify two different 
time-scales describing the system’s behavior: a fast time scale t where the regular oscillations occur and 
a slow time scale T = et which describes the slowly varying amplitude. Fig. [s)( A) depicts the variation 
of /?, which determines the oscillations frequency, with respect to the parameters r and B. Notice that 
/3 £ (0,1), i.e. /? oc 0(1). 
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Figure 4: Stability analysis of the spiral point L with respect to parameter B. (A,B) Stability 
analysis of L with respect to the functional tumor-associated vasculature B for the immune recruitment 
rate r equal to 0.54 and 0.60 days -1 , respectively. The labels S± and S 2 represent the bifurcation points, 
while the solid and dashed lines depict the stable and unstable solution branches, respectively. (C) Real 
part of the eigenvalues of the Jacobian estimated at the spiral point L with respect to B for r = 0.54 
days -1 . (D) Real part of the eigenvalues of the Jacobian estimated at the spiral point L with respect to 
B for r = 0.60 days -1 . 



0.55 0.65 0.75 0.85 

Recruitment rate r (days 1 ) 


Figure 5: Frequency dependence on model parameters B and r. Imaginary part of the eigenvalues 
of the Jacobian estimated at the spiral point L with respect to the immune recruitment rate r and the 
functional tumor-associated vasculature B. The labels S and U represent the regimes where the spiral 
point L is stable and unstable, respectively. 
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A plausible question relates to the possibility of controlling tumor growth and with the conditions 
under which this control can be achieved. The local stability analysis reveals that, when L is a spiral 
sink, the homoclinic orbit creates a trapping region where the tumor remains controlled. However, when 
L is a spiral source, the tumor will finally escape innate immune control, albeit the fact that it will stay 
to a certain region for a long period of time. Therefore, a potential strategy aiming at controlling the 
tumor growth would be to keep the tumor near the L fixed point. In the next sections, we show how an 
external modulation of the immune system dynamics could be designed to limit the uncontrolled tumor 
growth in the case that L is a spiral source. 


Mapping to a negatively damped harmonic oscillator via multi¬ 
scale analysis 


The original model given by Eqs. ®-(§ cannot be solved analytically and no control strategy can be 
easily designed. For this reason, we employ a multiscale approach to analyze the system’s behavior 
near the spiral point L. This approach, that efficiently describes the dynamics near a Hopf bifurcation 
point, reveals the inherent multiscale structure of the problem by capturing the regular oscillations and 
constructing an envelope of the slowly varying amplitude in deterministic 73 or stochastic nonlinear 
systems 74 76 . Thus, adopting such an approach we are able to map the initial system to a simplified 


negatively damped harmonic oscillator which allows to describe the self-sustained oscillations and the 
system’s energy gain 


77 


The system’s behavior near the spiral source L can be described by linearizing the Eqs. 0 ® around 
R = Rl and E = El. Then, we obtain a new system of equations approximating the evolution of the 
perturbations u and v given by: 


u = R — Rl , v = E — El, (3) 

for |rt| -C 1 and |t>| <C 1. Then, the linearized system takes the following form: 


d 

dt 


r of 



LdR 


L 

dF 1 1 

dE \ L 


ai 


L 

dG 1 
dE \ L - 


_«3 

a4 


( 4 ) 


where J|r represents the Jacobian at the fixed point L and F, G correspond to the right-hand side of 
Eqs. ©-©, respectively. The eigenvalues of the Jacobian are A = Tr(3\ L )/2 ± \J (Tr(J|i)) 2 — AD/2 = 
/r ± ip. Moreover, Tr(J|i) = 01+04 and D = 0404 — 0203 are the trace and determinant of the 
Jacobian matrix, respectively. Since e = /i and /3 oc 0(1), it follows that ai oc O(e), 04 oc O(e), 
whereas 02 oc 0(1) and 03 oc 0(1). Moreover, we observe that 0404 oc 0(e 2 ) which, consequently, 
results in fd = yj (ai + 04) 2 — 4(0404 — 0203 ) /2 ss yj—aia 3 . Thus, the eigenvalues can be approximated 
by A « A = e±zw, where u = yj— 0203 . The condition Kl holds for each case where L is a spiral source 
and 02 is always negative since it represents the death rate of tumor cells. Therefore, we can approximate 
the linearized system Q by the following one which incorporates the different order-terms: 


du 

dt 

dv 

dt 


eu + o 2 £, 
a 3 u + ev, 


( 5 ) 


where, u ~ u and v ~ v. The approximate linearized system © explicitly captures the different-order 
terms and allows to draw analogies with the field of Mechanics. More precisely, the system © represents a 
perturbed harmonic oscillator. The unperturbed problem (i.e. when e = 0) is a linear harmonic oscillator 
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with frequency w. The higher-order terms (or correction terms) insert small perturbations which result 
in a weakly negatively damped harmonic oscillator. 

The solution of the approximate system © can be expressed as: 

(j)=w(^) +B(T) {^n)’ (6 > 

where A[T) and B(T) contain the information regarding the slowly varying amplitude depending on T, 
where T = et for e <C 1. The multiscale assumption is that the functions A(T) and B(T ) evolve at the 
slow time scale T and are considered to be constant with respect to the oscillations with frequency w, 
evolving on the fast time scale t. Notice that the variables T = et oc O(e) and t oc 0(1) are considered to 
be independent. The constant term a 2 /u> in the expression of u has been used to simplify the upcoming 
calculations. 

In order to estimate the functions A(T) and B(T ), we assume that they follow evolution equations of 
the form: 


dA = f\dT, dB = f 2 dT. 


(7) 


Then, by taking the differentials of the system of equations in ([b]), we have that: 


du , du , . du 

d “= dt dt+ dA dA+ dB dB 


~ aivdt + (/, — sin(wt) + f 2 — cos(w()) dT, 


,, dv , dv , , 

dv= Si dt + aA dA + dB dB 


( 8 ) 


= a^udt + (/i cos(wt) — f -2 sin(wf)) dT. 

Notice that the system of equations in ([5]) and Q are equivalent. Therefore, by equating the terms which 
are of the same order, we obtain that: 


h 


02 

LO 


sin(o;t) + fi — cos (ut) = A(T)— sin(wf) + B{T)— cos {ojt), 

UJ LU UJ 

fi cos (ut) — f 2 sin(wt) = A(T ) cos (ujt) — B(T ) sin(wt). 


(9) 

( 10 ) 


To estimate the functions /i and f 2 , we project the Eqs. ([9|-(10l onto the fast dynamics. This is an 
important step of the calculations in order to isolate the amplitude of functions /i and f 2 75 76 . For 
instance, multiplying Eq. ([9]) by sin(wf) and integrating over a period of time [0, 2i r/w] , we have that: 



which results in /j = A. In a similar way, multiplying Eq. ([9]) by cos {uit) and following the same procedure, 
we find f 2 = B. Consequently, the solution of the approximate linearized system ([5]) is given by: 



= A 0 e T 


f sinM)\ 

cos (uit) J 


+ B 0 e T 


^ cos(wt)\ 
— sin (cot) J ’ 


( 11 ) 


where the positive parameters Aq and Bq are defined by the initial conditions of the original model 
©-©>■ In addition, Eq. (|TT|) can be further simplified to take the following form: 
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fMe T sin(wt + <j))\ 

\v J yNe T sin(wf + 8)) ' 

(12) 

where M cos(<f>) = A 0 a 2 /co, M sin(</>) = B 0 a 2 /uj, N sm(S) = A 0 and N cos(5) = 
functions u and v are orthogonal since: 

—Bo. We show that the 

cos(<5 — (j>) = cos (8) cos((/)) + sin(<5) sin(^) 

(13) 

A) «2 Ao Bq a 2 

N M uj + N M oj _ 

(14) 

r , 7T 

^8 = 0+ 

(15) 


Then, the system of the approximate solutions becomes: 


/ Me T sin(o>t + <j>) 
\Ne T cos (ut + (j>) 


(16) 


The advantage of using the proposed approach is that the solutions of the approximate model (161 


depend directly on experimentally accessible parameters and the initial conditions. We focus now on the 
stability properties of the system (|5j) with respect to the time evolution of the variables u and v. To that 
end, we construct a Lyapounov functional as: 


(17) 


V (u, v) = ^ (a 3 u 2 - a 2 v 2 ) , 

where V is always positive definite since a 2 < 0, and the time-derivative of V is given by: 

d V „ du , dv 

~di = a3U M ~ U2V dt 

= e (a 3 u 2 — a 2 v 2 ) + {aza 2 uv — asa 2 uv ) 

= e ( a 3 u 2 - a 2 v 2 ) > 0. 

Since dV/dt > 0 the system gains energy. The average energy (V) over a period is estimated by using 


(18) 


the form of the approximate solutions (121 and considering the multiscale assumption: 

r 2n/uj i i 

(V) = Vdt = - ( a 3 u 2 — a 2 D 2 ) dt « e 2T - (a 3 M 2 — a 2 iV 2 ) . 

Jo 2 4 

The term | (a 3 M 2 — a 2 lV 2 ) is constant and depends on time-invariant parameter values and the initial 
conditions. Consequently, the average energy gain rate over a period is equal to: 


(19) 


d(V) 

dt 


= MV). 


( 20 ) 


The previous results suggest that a therapeutic strategy, which influences the effector cell dynamics, 
should be designed in a way to “pump out energy” with an average rate which is greater than the system’s 
gain energy. In particular, the per period gain rate is equal to 2e according to relation (191. Thus, if we 
introduce an external immune-modulatory term h(v) to intervene in the tumor-effector cell interactions, 
the system of equations (J5| becomes: 


du 

dt 

dv 

dt 


eu + a 2 v , 
a 3 u + ev + h(v). 


( 21 ) 
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According to Eq. ( |l8| ), the rate dV/dt is at most of the order of 0(e). Therefore, the function h(v) should 
be of the order of 0(e), for instance: 


h(v) = cz(t). 


( 22 ) 


In order to calculate the function z(t), we require that the energy 14 of the system (21) should have 
a negative rate, that is: 


dVh 


du 


Ilf = a3U dt 


a 2 v 


, dv 
dt 


= e (a^u 2 — a 2 v 2 — a 2 vz(t )) + {a^,a 2 uv — a^a 2 uv) 
= e {a^u 2 — a 2 v 2 — a 2 vz(tj) < 0. 


(23) 


The negative energy rate means that a trajectory will flow “downhill” towards a stable fixed point. In 
the limit dVh/dt = 0, we find an expression for the zero-rate energy function z(t) in terms of the solutions 
u and v, given by: 


z{t) = 


a^u 2 — a 2 i’ 2 
a 2 v 


(24) 


Substituting the solutions of u and v from Eq. (16), we have an analytical expression of the zero-rate 
energy function z(t), that is: 


z(t) = 


U3 114“ rp 

a 2 N 6 


tan (cat + </>) sin(wl + <j>) 


Me T cos (cot + (f>). 


(25) 


The function z(t) has singularities at the points where the function cos (cat + <p) is equal to zero, i.e. 
tsing = ( K7r + f — <^) /w, k = 0,1,2,.... However, the effect of z(t), i.e. the integral of the function within 
a small time interval [ti : t 2 ] C [0,27r/w], is finite. 

Fig. [6] illustrates the effect of z(t) by integrating between t = 0 and t = tt/w, for values of the model 
parameters r and B where L is unstable. In each case, the system was initialized near the spiral point L. 



Recruitment rate r (days 1 ) 


Figure 6: Dependence of the zero-rate energy function z(t) effect on model parameters B and 

r. Dependence of the effect of z(t) on the immune recruitment rate r and functional tumor-associated 
vasculature B in the unstable regime of L denoted by U. The label S stands for the regime where the 
spiral point is stable and tumor control is reached. 
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Using the function z(t) is not feasible for practical purposes. However, we can design meaningful 
therapy functions z ext {t) that induce the same or greater effects in a certain time interval [t \, t 2 ] as the 
function z(t), that is: 


z(t) dt < / z ext (t)dt. 


(26) 


This will result in “pumping out energy” at a rate greater than the system’s energy gain. In the 
next section, we estimate the function z(t) which represents the zero-rate energy scenario, as well as we 
present numerical simulation results for specific values of the model parameters. Furthermore, based on 
the behavior of the function z(t), we suggest an efficient external immune-modulatory term z ex t(t ) to 
fulfill the relation (26) which results in long-term tumor control. 


Results: theory-driven therapeutic design 


In this section, we design an immuno-tlierapeutic proposal derived from our model analysis. We first 
compare the approximate solutions with those obtained from the original model 0-0- Then, we 
show how the defined energy function in Eq. 0 for the approximate linearized system can be used to 
describe the system’s energy gain. Finally, we design an external immuno-modulatory strategy following 


the behavior of the zero-rate energy function given in Eq. (251. 


To illustrate the simulation results and without loss of generality, we select the following parameter 
values: r = 0.6 days -1 and B = 0.8, which result in e = 5.85 • 10 -4 and w = ^/—a 2 a 3 = 0.336 rad/days. 
Notice that similar results can be obtained for each parameter set where L is a spiral source, since eCl 
always holds. For numerical integration, we consider a d th order Runge-Kutta method. The time step 
was set equal to dt = 0.001 days which is sufficiently small. 


Validity of approximate solution 

At this point, we need to validate the performance of the approximate solutions. Fig. [7)( A) depicts the 
evolution of the approximate perturbation u compared to the perturbation u = R — Rl of the initial 
model given by Eqs. 0-0- referred in what follows as original model perturbation. The approximate 
solutions u were estimated by performing a numerical integration of the system ([5|, while u = R — Rl 
by numerically integrating the original model 0-©. Both systems of equations were initialized close 
to the spiral point L and the simulation time was set t S i m = 1000 days. A quick glance at Fig. 0A) 
reveals that the solution derived by the approximate model ([5]) is very close to the original one, which 
demonstrates the validity of the approach. Moreover, Fig. [7](A) shows a comparison of the original and 
approximate model solutions by zooming in the time interval between days 400 and 425. In this interval, 
the maximum error between such solutions was found to be of the order of 10 -3 . 
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Figure 7: Validation of approximate solutions. (A) Evolution of the approximate perturbation u 
compared to the original model perturbation u = R — Rl■ (B) Evolution of the energy function of 
the approximate linearized model defined in Eq. ( |17| against the energy estimated by using the original 
model 0-0- (C) The zero-rate energy function in Eq. ( |24| ) estimated by using the approximate solution 
of the system ©. (D) The immuno-modulatory function z ex t(t ) for a = —0.6 and 7 = 100 compared 
with the zero-rate energy function z(t). 


Fig. [tJ[B) compares the evolution of the approximate energy function defined in Eq. ( [17) against the 
energy of the original model As in the previous case, the energy function of the approximation is 

in a good agreement with that obtained from the original model. Therefore, we expect that an external 
immune-modulatory strategy, which results in a negative energy rate, can be used with the same efficacy 
to the original model. 


Construction of the therapy term 

Fig. [7)C) depicts the zero-rate energy function in Eq. (24) estimated by using the approximate solutions 
u and v from the system ( 21 ). As we observe, this function diverges at the singularity points, which 


means that it cannot be used directly to induce tumor control. 

It should be noted that, the zero-rate energy function changes sign according to the evolution of the 
variable u = R — Rl , he. the function follows the monotonicity of tumor evolution. Hence, an external 


immuno-modulatory therapy z ex t ( t ) that satisfies the condition (261 should follow the evolution of the 


tumor by increasing the recruitment rate of effector cells when the tumor radius increases and decreasing 
the immune recruitment rate when the tumor radius decreases. A simple approximation of the therapy 
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function would be a step function having a constant value and the same sign to the zero-rate energy 
function at the same time interval: 


z e xt(t) = a tanh( 7 #), (27) 

where 7 > 1 expresses the sigmoidal steepness, and the selection of a should fulfill the relation ( [26] ). 
Fig.^D) depicts the immuno-modulatory function z ex t(t) for a = — 0.6 and 7 = 100 , compared with the 
zero-rate energy function z(t) estimated. 

Performance of the suggested therapy 

The application of the external therapy in Eq. ( |27| ) implies the addition of a new term in Eqs. 0 -©. 
that is: 


dR 

dt 


-(XmB — A a)R + Am(1 - B)L d 




\tanh(R/Lr>) R 


- — cERf(R, B), 


dE 


R? 


dt r K + R 3 


E — d\ER? f(R,a) — doE + a + eatanh( 7 u(t)), 


(28) 

(29) 


where v(t) = E(t) — E^. 

Fig. | 8 ]A,B) respectively shows the evolution of the tumor radius and effector cells by numerically 
integrating Eqs. (28)-(291 which include the proposed external immuno-modulatory function (27). The 
system is initialized near the spiral point L and parameter values of the external function were selected 
equal to a = —0.6 and 7 = 100. In doing so, we obtain that the system reaches a stable steady state, 
i.e. the tumor remains controlled. Fig. | 8 [C) illustrates how the energy of the system of Eqs. (28)-(29) 
evolves. The energy decreases with time and becomes zero as the system reaches the steady state. It 
is worth pointing out that a was selected not only to satisfy the relation (26), but also to result in the 
system’s fast convergence to a steady state. 
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Figure 8 : Performance of the proposed therapy term. Evolution of the tumor radius (A) and 
concentration of effector cells (B) by considering the external immuno-modulatory function (27) with 
parameters a = —0.6 and 7 = 100. (C) Energy of the system of Eqs. (28)-(29) with a = —0.6 and 
7 = 100. Evolution of the tumor radius (D) and concentration of effector cells (E) by initializing the 
system away from the spiral point L and considering the external immuno-modulatory function with 
a = —25 and 7 = 100. (F) Energy of the system of Eqs. (28)-(291, with a = —25 and 7 = 100, initialized 
away from the spiral point L. 


Fig. [ 8 j^D,E) respectively shows the evolution of the tumor radius and concentration of effector cells 
by numerically integrating Eqs. (28)-(29) when the system is initialized away from the spiral point L. 
In this case, the approximate solution is expected to deviate from the solutions of the original model. 
Fig.^F) represents the temporal evolution of the corresponding energy. The parameter a was selected 
to be equal to —25 to fulfil the relation (26), as well as to provide a fast convergence to a steady state. 
Interestingly enough, in this case the system also reaches a steady state, even though the approximate 
system is not accurate enough. Consequently, the proposed external immuno-modulatory function is 
shown to be adequate in controlling tumor growth, even when the system is initialized far to the spiral 
point L. 


Discussion 


In this article, we investigate the therapeutic potential of immunomodulatory interventions against tumor 


growth. To that end, we consider a model introduced in 53 that describes the dynamic interplay between 


vascularized tumor growth and effector cell responses. Our goal is to identify an external modification 
of effector cell dynamics that allows for controlling tumor growth. With the help of bifurcation analysis, 
we identified a unstable oscillatory regime that induces tumor evasion from immuno-surveillance. The 
characteristic feature of this regime is the oscillations occur at a faster time scale than the amplitude 
dynamics. Exploiting this time scale separation and via temporal multiscale analysis, we map our model 
onto a weakly negatively damped harmonic oscillator. This approximation allowed us to identify an 
analytical expression for the additive control term to the effector cell dynamics. This term acts as an 
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external immunomodulatory therapy where the numerical simulations evidence its efficacy in controlling 
tumor growth. 

The crucial question concerns the translational potential of our theory-driven therapeutic proposal into 
clinical practice. Our suggested immunomodulatory strategy is relevant to small enough, non-invasive 
tumors that are initially controlled by the immune system but eventually evade immuno-surveillance. 
The latter occurs when tumors exceed a critical size where immune responses are unable to confer any 
control. As stated above, our model suggests that such tumor evasion may take place in the form of 
oscillations of slowly increasing amplitude. The proposed therapeutic strategy is based on the syn¬ 
chronization of immuno-stimulating and -suppressive phases with tumor growth dynamics. Although 
immuno-stimulating therapies seem to be expected and plausible |;6,[8], immuno-suppression sounds 
counter-intuitive and dangerous. However, the latter occurs in clinical practice during chemotherapeutic 
interventions 78 . Therefore, a potential realization of our proposed strategy could be mediated by a com¬ 
bination of state-of-the-art immunotherapies [l|l4l|8,49 and chemotherapeutic modules. The latter not 


only will play the role of immuno-suppressor, but also will slow down tumor growth dynamics. Needless 
to state that such a therapeutic suggestion requires experimental validation and further investigation. 

Finally, we conclude by pointing out the limitations of the present work. Paramount among them, 
introducing vascularization dynamics should be not only a natural, but also an insightful extension of 
the proposed model. At this state tumor vascularization B is considered as a constant parameter in 
time. Thus, making vascularization dynamic would make sense to model further hypoxic effects, such as 
necrosis or hypoxic-induced invasion. Moreover, the immune system is much more complicated than its 
description in the current model, involving much more cell types and interactions. In particular, immune 
system is often regarded as acting to suppress tumor growth, however it is now clear that it can be both 
stimulatory and inhibitory, as in the case of tumor-associated macrophages 79 80 . Moreover, our model 
needs be extended for invasive tumors. Nevertheless, mathematical modeling can help to improve our 
understanding of the interplay between the several competing factors that have complex implications for 
cancer therapy. We hope that the model presented here could provide a step towards this direction. 


Acknowledgments 

A. I. Reppas, J. C. L. Alfonso and H. Hatzikirou gratefully acknowledge the support of the German 
Excellence Initiative via the Cluster of Excellence EXC 1056 Center for Advancing Electronics Dresden 
(cfAED) at the Technische Universitat Dresden. J. C. L. Alfonso and H. Hatzikirou also acknowledge 
the German Federal Ministry of Education and Research (BMBF) for the eMED project SYSIMIT 
(01ZX1308D). 

References 

[1] O. Finn, “Immuno-oncology: understanding the function and dysfunction of the immune system in 
cancer,” Ann. Oncol., vol. 23, no. 8, pp. 6—9, 2012. 

[2] G. Dranoff, “Cytokines in cancer pathogenesis and cancer therapy,” Nat. Rev. Cancer, vol. 4, no. 1, 
pp. 11 22, 2004. 

[3] T. Stewart and S. Abrams, “How tumours escape mass destruction,” Oncogene, vol. 27, no. 45, 
pp. 5894-903, 2008. 

[4] W. Zou, “Immunosuppressive networks in the tumour environment and their therapeutic relevance,” 
Nat. Rev. Cancer, vol. 5, no. 4, pp. 263-74, 2005. 








18 


[5] G. Rabinovich, D. Gabrilovich, and E. Sotomayor, “Immunosuppressive strategies that are mediated 
by tumor cells,” Annu. Rev. Immunol, vol. 25, pp. 267-96, 2007. 

[61 C. Ghirelli and T. Hagemann, “Targeting immunosuppression for cancer therapy,” J. Clin. Invest., 
vol. 123, no. 6, pp. 2355-7, 2013. 

[7] M. Cheever, J. Seldom, L. Weiner, H. Lyerly, M. Disis, A. Greenwood, O. Grad, and W. Nelson, 
“Translational research working group developmental pathway for immune response modifiers,” Clin. 
Cancer Res., vol. 14, no. 18, pp. 5692-9, 2008. 

[8] S. Rosenberg, “Decade in review-cancer immunotherapy: entering the mainstream of cancer treat¬ 
ment,” Nat. Rev. Clin. Oncol, vol. 11, no. 11, pp. 630-2, 2014. 

[9] O. Finn, “Cancer immunology,” N. Engl. J. Med,., vol. 358, no. 25, pp. 2704-15, 2008. 

[10] D. Ambrosi and F. Mollica, “On the mechanics of a growing tumor,” International journal of engi¬ 
neering science, vol. 40, no. 12, pp. 1297 1316, 2002. 

[11] H. Byrne and L. Preziosi, “Modelling solid tumour growth using the theory of mixtures,” Mathe¬ 
matical Medicine and Biology, vol. 20, no. 4, pp. 341 -366, 2003. 

[12] R. K. Jain, J. D. Martin, and T. Stylianopoulos, “The role of mechanical forces in tumor growth 
and therapy,” Annual review of biomedical engineering, vol. 16, p. 321, 2014. 

[13] A. Ramirez-Torres, R. Rodriguez-Ramos, J. Merodio, J. Bravo-Castillero, R. Guinovart-Diaz, and 
J. Alfonso, “Action of body forces in tumor growth,” International Journal of Engineering Science, 
vol. 89, pp. 18-34, 2015. 

[14] A. Ramirez-Torres, R. Rodriguez-Ramos, J. Merodio, J. Bravo-Castillero, R. Guinovart-Diaz, and 
J. Alfonso, “Mathematical modeling of anisotropic avascular tumor growth,” Mechanics Research 
Communications, vol. 69, pp. 8-14, 2015. 

[15] J. Arciero, T. Jackson, and D. Kirschner, “A mathematical model of tumor-immune evasion and 
sirna treatment,” Discrete Continuous Dyn Syst Ser B , vol. 4, pp. 39-58, Feb. 2004. 

[16] N. Bellomo, B. Firmani, and L. Guerri, “Bifurcation analysis for a nonlinear system of integro- 
differential equations modelling tumor-immune cells competition,” Appl Math Lett, vol. 12, pp. 39- 
44, Mar. 1999. 

[17] N. Bellomo and N. Preziosi, “Modelling and mathematical problems related to tumor evolution and 
its interaction with the immune system,” Math. Comput. Modelling, vol. 32, pp. 413-452, Aug. 2000. 

[18] M. Galach, “Dynamics of the tumor-immune system competition - the effect of time delay,” Int J 
Appl Math Comput Sci, vol. 13, no. 3, pp. 395-406, 2003. 

[19] A. Matzavinos, M. Chaplain, and V. Kuznetsov, “Mathematical modelling of the spatio-temporal 
response of cytotoxic t-lymphocytes to a solid tumour,” Math Med Biol, vol. 21, pp. 1 34, Mar. 2004. 

[20] A. d’Onofrio, “A general framework for modeling tumor-immune system competition and im¬ 
munotherapy: Mathematical analysis and biomedical inferences,” Physica D: Nonlinear Phenomena, 
vol. 208, pp. 220-235, Sept. 2005. 

[21] D. Mallet and L. De Pillis, “A cellular automata model of tumor-immune system interactions,” J 
Theor Biol, vol. 239, pp. 334-350, Apr. 2006. 



19 


[22] M. Al-Tameemi, M. Chaplain, and A. d’Onofrio, “Evasion of tumours from the control of the immune 
system: consequences of brief encounters,” Biol Direct, vol. 7, p. 31, Sept. 2012. 

[23] M. Robertson-Tessi, A. El-Kareh, and A. Goriely, “A mathematical model of tumor-immune inter¬ 
actions,” J Theor Biol, vol. 294, pp. 56-73, Feb. 2012. 

[24] F. Rihan, D. Abdel Rahman, S. Lakshmanan, and A. Alkhajeli, “A time delay model of tumour- 
immune system interactions: Global dynamics, parameter estimation, sensitivity analysis,” Appl 
Math Comput, vol. 232, pp. 606-623, Apr. 2014. 

[25] C. Koebel, W. Vermi, J. Swann, N. Zerafa, S. Rodig, L. Old, M. Smyth, and R. Schreiber, “Adaptive 
immunity maintains occult cancer in an equilibrium state,” Nature, vol. 450, pp. 903-907, 2007. 

[26] A. d’Onofrio, “Tumor evasion from immune control: Strategies of a miss to become a mass,” Chaos 
Soliton. Fract., vol. 31, no. 2, pp. 261-268, 2007. 

[27] G. Caravagna, A. d’Onofrio, P. Milazzo, and R. Barbuti, “Tumour suppression by immune system 
through stochastic oscillations,” J Theor Biol, vol. 265, pp. 336-345, Aug. 2010. 

[281 D. Pardoll, “Does the immune system see tumors as foreign or self?,” Annu. Rev. Immunol., vol. 21, 
pp. 807-39, 2003. 

[29] G. Dunn, L. Old, and R. Schreiber, “The three es of cancer immunoediting,” Annu. Rev. Immunol., 
vol. 22, pp. 329-60, 2004. 

[301 J. Stark, C. Chan, and A. George, “Oscillations in the immune system,” Immunol Rev, vol. 216, 
pp. 213-231, Apr. 2007. 

[31] B. Coventry, M. Ashdown, M. Quinn, S. Markovic, S. Yatomi-Clarke, and A. Robinson, “Crp iden¬ 
tifies homeostatic immune oscillations in cancer patients: a potential treatment targeting tool?,” J. 
Transl. Med,., vol. 7, p. 102, 2009. 

[32] V. Kuznetsov, I. Makalkin, M. Taylor, and A. Perelson, “Nonlinear dynamics of immunogenic tumors: 
Parameter estimation and global bifurcation analysis,” Bull Math Biol, vol. 56, pp. 295-321, Mar. 
1994. 

[33] D. Kirschner and J. Panetta, “Modeling immunotherapy of the tumor-immune interaction,” J Math 
Biol , vol. 37, pp. 235-52, Sept. 1998. 

[34] A. d’Onofrio, “Tumour-immune system interaction: Modeling the tumour-stimulated proliferation 
of effectors and immunotherapy,” Math Models Methods Appl Sci, vol. 16, pp. 1375-1401, Aug. 2006. 

[35] L. de Pillis, W. Gu, and A. Radunskaya, “Mixed immunotherapy and chemotherapy of tumors: 
modeling, applications and biological interpretations,” J Theor Biol, vol. 238, pp. 841-862, Feb. 
2006. 

[36] A. d’Onofrio, F. Gatti, P. Cerrai, and L. Freschi, “Delay-induced oscillatory dynamics of tumour- 
immune system interaction,” Math Comput Model, vol. 51, pp. 572-591, Mar. 2010. 

[37] R. Araujo and D. McElwain, “A history of the study of solid tumour growth: the contribution of 
mathematical modelling,” Bull Math Biol, vol. 66, pp. 1039-1091, Sept. 2004. 

[38] J. Nagy, “The ecology and evolutionary biology of cancer: a review of mathematical models of 
necrosis and tumor cell diversity,” Math Biosci Eng, vol. 2, pp. 381-418, Apr. 2005. 



20 


[39] H. Byrne, T. Alarcon, M. Owen, S. Webb, and P. Maini, “Modelling aspects of cancer dynamics: a 
review,” Philos Trans A Math Phys Eng Sci , vol. 364, pp. 1563-1578, June 2006. 

[40] T. Roose, S. Chapman, and P. Maini, “Mathematical models of avascular tumor growth,” SIAM 
Rev, vol. 49, pp. 179-208, May 2007. 

[41] M. Martins, S. Ferreira Jr., and M. Vilela, “Multiscale models for the growth of avascular tumors,” 
Phys Life Rev, vol. 4, pp. 128-156, June 2007. 

[42] N. Bellomo and M. Delitala, “From the mathematical kinetic, and stochastic game theory to mod¬ 
elling mutations, onset, progression and immune competition of cancer cells,” Phys Life Rev, vol. 5, 
pp. 183-206, Dec. 2008. 

[43] M. Chaplain, “Modelling aspects of cancer growth: Insight from mathematical and numerical anal¬ 
ysis and computational simulation,” in Multiscale Problems in the Life Sciences (V. Capasso and 
M. Lachowicz, eds.), vol. 1940 of Lecture Notes in Mathematics, pp. 147-200, Springer Berlin Hei¬ 
delberg, 2008. 

[44] R. Eftimie, J. Bramson, and D. Earn, “Interactions between the immune system and cancer: a brief 
review of non-spatial mathematical models,” Bull Math Biol, vol. 73, pp. 2-32, Jan. 2011. 

[45] K. Wilkie, “A review of mathematical models of cancer-immune interactions in the context of tumor 
dormancy,” Adv Exp Med Biol, vol. 734, pp. 201-234, Aug. 2013. 

[46] E. Ruoslahti, “Specialization of tumour vasculature,” Nat Rev Cancer, vol. 2, pp. 83-90, Feb. 2002. 

[47] N. Ferrara and R. Kerbel, “Angiogenesis as a therapeutic target,” Nature, vol. 438, pp. 967-974, 
Dec. 2005. 

[48] R. Farnsworth, M. Lackmann, M. Achen, and S. Stacker, “Vascular remodeling in cancer,” Oncogene, 
vol. 33, pp. 3496-3505, Aug. 2014. 

[49] W. Fridman, F. Pages, C. Sautes-Fridman, and J. Galon, “The immune contexture in human tu¬ 
mours: impact on clinical outcome,” Nat Rev Cancer, vol. 12, pp. 298-306, Mar. 2012. 

[50] M. Junttila and F. de Sauvage, “Influence of tumour micro-environment heterogeneity on therapeutic 
response,” Nature, vol. 501, pp. 346-354, Sept. 2013. 

[51] H. Greenspan, “On the growth and stability of cell cultures and solid tumors,” J Theor Biol, vol. 56, 
pp. 229-242, Jan. 1976. 

[52] V. Cristini, J. Lowengrub, and Q. Nie, “Nonlinear simulation of tumor growth,” J Math Biol, vol. 46, 
pp. 191-224, Mar. 2003. 

[53] H. Hatzikirou, J.C.L. Alfonso, S. Miihle, C. Stern, S. Weiss, and M. Meyer-Hermann, “Cancer thera¬ 
peutic potential of combinatorial innnuno- and vaso-modulatory interventions,” arXiv:1505.05670vl, 
2015. 

[54] H. Hatzikirou, A. Chauviere, J. Lowengrub, J. De Groot, and V. Cristini, “Effect of vascularization 
on glioma tumor growth,” in Modeling Tumor Vasculature (T. L. Jackson, ed.), pp. 237-259, Springer 
New York, 2012. 

[55] L. de Pillis, A. Radunskaya, and C. Wiseman, “A validated mathematical model of cell-mediated 
immune response to tumor growth,” Cancer Res, vol. 65, pp. 7950-8, Sept. 2005. 


21 


[56] Y. Huang, S. Goel, D. Duda, D. Fukumura, and R. Jain, “Vascular normalization as an emerging 
strategy to enhance cancer immunotherapy,” Cancer Res , vol. 73, pp. 2943-2948, May 2013. 

[57] A. Herman, V. Savage, and G. West, “A quantitative theory of solid tumor growth, metabolic rate 
and vascularization,” PLoS ONE, vol. 6, no. 9, p. e22973, 2011. 

[58] A. Bru, S. Albertos, J. Subiza, J. Garcia-Asenjo, and I. Bru, “The universal dynamics of tumor 
growth,” Biophys J, vol. 85, pp. 2948-61, Nov. 2003. 

[59] A. Bru, S. Albertos, J. Garcia-Asenjo, and I. Bru, “Pinning of tumoral growth by enhancement of 
the immune response,” Phys Rev Lett , vol. 92, pp. 1-4, June 2004. 

[60] K. Alessandri, B. Sarangi, V. Gurchenkov, B. Sinha, T. Kiefiling, L. Fetler, F. Rico, S. Scheming, 
C. Lamaze, A. Simon, S. Geraldo, D. Vignjevic, H. Domejean, L. Rolland, A. Funfak, J. Bibette, 
N. Bremond, and P. Nassoy, “Cellular capsules as a tool for multicellular spheroid production and 
for investigating the mechanics of tumor progression in vitro,” Proc Natl Acad Sci USA, vol. 110, 
pp. 14843-14848, Aug. 2013. 

[61] M. Delarue, F. Montel, O. Caen, J. Elgeti, J. Siaugue, D. Vignjevic, J. Prost, J. Joanny, and 
G. Cappello, “Mechanical control of cell flow in multicellular spheroids,” Phys Rev Lett, vol. 110, 
p. 138103, Mar. 2013. 

[62] H. Acker, Spheroids in cancer research: methods and perspectives. Springer-Verlag Berlin; New York, 
1984. 

[63] H. Frieboes, X. Zheng, C. Sun, B. Tromberg, R. Gatenby, and V. Cristini, “An integrated computa¬ 
tional/experimental model of tumor invasion,” Cancer Res, vol. 66, pp. 1597 1604, Feb. 2006. 

[64] V. Cristini, H. Frieboes, X. Li, J. Lowengrub, P. Macklin, S. Sanga, S. Wise, and X. Zheng, “Non¬ 
linear modeling and simulation of tumor growth,” in Selected topics in cancer modeling: Genesis, 
evolution, immune competition, and therapy. Modelling and Simulation in Science, Engineering, and 
Technology (N. Bellomo, M. Chaplain, and E. de Angelis, eds.), ch. 6, pp. 113-82, Boston, MA USA: 
Birkhauser, 2008. 

[65] B. Su, W. Zhou, K. Dorman, and D. Jones, “Mathematical modelling of immune response in tissues,” 
Comput Math Methods Med, vol. 10, no. 1, pp. 9-38, 2009. 

[66] A. d’Onofrio, U. Ledzewicz, and H. Schattler, “On the dynamics of tumor-immune system in¬ 
teractions and combined chemo- and immunotherapy,” in New Challenges for Cancer Systems 
Biomedicine (A. d’Onofrio, P. Cerrai, and A. Gandolfi, eds.), SIMAI Springer Series, pp. 249-266, 
Springer Milan, 2012. 

[67] C. Kelley, Iterative methods for optimization, vol. 18. SIAM, 1999. 

[68] H. B. Keller, “Numerical solution of bifurcation and nonlinear eigenvalue problems,” Applications of 
bifurcation theory, pp. 359-384, 1977. 

[69] A. Nayfeh and B. Balachandran, Applied nonlinear dynamics: analytical, computational and exper¬ 
imental methods. John Wiley and Sons, 2008. 

[70] S. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry and 
engineering. Perseus Books Group, 2001. 

[71] W. E. Boyce, R. C. DiPrima, and C. W. Haines, Elementary differential equations and boundary 
value problems, vol. 9. Wiley New York, 1992. 



22 


[72] G. Iooss and D. D. Joseph, Elementary stability and bifurcation theory. Springer Science & Business 
Media, 2012. 

[73] J. Cole and J. Kevorkian, “Multiple scale and singular perturbation methods,” Appl. Math. Sci, 
vol. 114, 1996. 

[74] M. Klosek and R. Kuske, “Multiscale analysis of stochastic delay differential equations,” Multiscale 
Model Simul., vol. 3, no. 3, pp. 706-729, 2005. 

[75] R. Kuske, “Multi-scale analysis of noise-sensitivity near a bifurcation,” in IUTAM Symposium on 
Nonlinear Stochastic Dynamics, pp. 147-156, Springer Netherlands, 2003. 

[76] R. Kuske, L. Gordillo, and P. Greenwood, “Sustained oscillations via coherence resonance in sir,” J 
Theor Biol, vol. 245, no. 3, pp. 459-469, 2007. 

[77] A. Jenkins, “Self-oscillation,” Physics Reports, vol. 525, no. 2, pp. 167 -222, 2013. 

[78] L. Zitvogel, L. Apetoh, F. Ghiringhelli, and G. Kroemer, “Immunological aspects of cancer 
chemotherapy,” Nat Rev Immunol, vol. 8, no. 1, pp. 59-73, 2008. 

[79] K. de Visser, A. Eichten, and L. Coussens, “Paradoxical roles of the immune system during cancer 
development,” Nat Rev Cancer, vol. 6, pp. 24-37, Jan. 2006. 

[80] A. Mantovani, P. Romero, A. Palucka, and F. Marincola, “Tumour immunity: effector response to 
tumour and role of the microenvironment,” Lancet, vol. 371, pp. 771-783, Mar. 2008. 



